Hierarchical model

Rats experiment

(BDA3, p. 102)

# part of the code below is copied from 
# Aki Vehtari <Aki.Vehtari@aalto.fi>
# Markus Paasiniemi <Markus.Paasiniemi@aalto.fi>

library(ggplot2)
library(gridExtra)
library(tidyr)

y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
       2,1,5,2,5,3,2,7,7,3,3,2,9,10,4,4,4,4,4,4,4,10,4,4,4,5,11,12,
       5,5,6,5,6,6,6,6,16,15,15,9,4)
n <- c(20,20,20,20,20,20,20,19,19,19,19,18,18,17,20,20,20,20,19,19,18,18,25,24,
       23,20,20,20,20,20,20,10,49,19,46,27,17,49,47,20,20,13,48,50,20,20,20,20,
       20,20,20,48,19,19,19,22,46,49,20,20,23,19,22,20,20,20,52,46,47,24,14)

x <- seq(0.0001, 0.9999, length.out = 1000)
summary(y/n)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.05132 0.11111 0.13812 0.21053 0.37500
hist(y/n)

Fit beta-binomial with flat prior separately

Each experiment has a different theta

# helperf function to evaluate density over observations
bdens <- function(n, y, x)
  dbeta(x, y+1, n-y+1)

df_sep <- mapply(bdens, n, y, MoreArgs = list(x = x)) %>%
  as.data.frame() %>% cbind(x) %>% gather(ind, p, -x)

labs1 <- paste('posterior of', c('theta_j', 'theta_71'))

# plot the separate model
plot_sep <- ggplot(data = df_sep) +
  geom_line(aes(x = x, y = p, color = (ind=='V71'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Separate model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
# The last one is for emphasize colored red
plot_sep

Fit pooled Beta-Binomial model (common theta)

df_pool <- data.frame(x = x, p = dbeta(x, sum(y)+1, sum(n)-sum(y)+1))

plot_pool <- ggplot(data = df_pool) +
  geom_line(aes(x = x, y = p, color = '1')) +
  labs(x = expression(theta), y = '', title = 'Pooled model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = 'red', labels = 'Posterior of common theta') +
  theme(legend.background = element_blank(), legend.position = c(0.7,0.9))

# Plot both separate and pooled model
grid.arrange(plot_sep, plot_pool)

Hierarchical model (Beta-Binomial)

A <- seq(0.1, 10, length.out = 200) ## alpha
B <- seq(3, 60, length.out = 200) ## beta
# make vectors that contain all pairwise combinations of A and B
cA <- rep(A, each = length(B))
cB <- rep(B, length(A))
# Use logarithms for numerical accuracy!
lpfun <- function(a, b, y, n) 
  sum(lgamma(a+b)-lgamma(a)-lgamma(b)+lgamma(a+y)+lgamma(b+n-y)-lgamma(a+b+n))

lp <- mapply(lpfun, cA, cB, MoreArgs = list(y, n))
df_marg <- data.frame(x = cA, y = cB, p = exp(lp - max(lp)))

# Subtract maximum value to avoid over/underflow in exponentation
title1 <- 'Contour of likelihood for alpha beta'
# create a plot of the marginal posterior density
postdensityalphabeta = ggplot(data = df_marg, aes(x = x, y = y)) +
  geom_raster(aes(fill = p, alpha = p), interpolate = T) +
  geom_contour(aes(z = p), colour = 'black', size = 0.2) +
  coord_cartesian(xlim = range(cA), ylim = range(cB)) +
  labs(x = 'alpha', y = 'beta', title = title1) +
  scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
  scale_alpha(range = c(0, 1), guide = F)
postdensityalphabeta

First visualize(calculate) the marginal posterior distribution of alpha and beta after integrating out theta’s

# Compute the marginal posterior of alpha and beta in hierarchical model
# Use grid
A <- seq(0.5, 6, length.out = 100) ## alpha
B <- seq(3, 33, length.out = 100) ## beta
# make vectors that contain all pairwise combinations of A and B
cA <- rep(A, each = length(B))
cB <- rep(B, length(A))

# Use logarithms for numerical accuracy!
lpfun <- function(a, b, y, n) log(a+b)*(-5/2) +
  sum(lgamma(a+b)-lgamma(a)-lgamma(b)+lgamma(a+y)+lgamma(b+n-y)-lgamma(a+b+n))

lp <- mapply(lpfun, cA, cB, MoreArgs = list(y, n))
df_marg <- data.frame(x = cA, y = cB, p = exp(lp - max(lp)))

# Subtract maximum value to avoid over/underflow in exponentation
title1 <- 'The marginal posterior of alpha and beta in hierarchical model'
# create a plot of the marginal posterior density
postdensityalphabeta = ggplot(data = df_marg, aes(x = x, y = y)) +
  geom_raster(aes(fill = p, alpha = p), interpolate = T) +
  geom_contour(aes(z = p), colour = 'black', size = 0.2) +
  coord_cartesian(xlim = c(1,5), ylim = c(4, 26)) +
  labs(x = 'alpha', y = 'beta', title = title1) +
  scale_fill_gradient(low = 'yellow', high = 'red', guide = F) +
  scale_alpha(range = c(0, 1), guide = F)
postdensityalphabeta

obtain random samples from p(alpha,beta|data)

(samp_A, samp_B) contains 100 samples of alpha and beta

# df_marg: first column (value of alpha)
#          second column (value of beta)
#          third column (value of p(alpha, beta | data))
# sample from the grid (with replacement)
nsamp <- 100
samp_indices <- sample(length(df_marg$p), size = nsamp,
                       replace = T, prob = df_marg$p/sum(df_marg$p))

samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]

visualize the samples

nsamp2 <- 2000
samp_indices2 <- sample(length(df_marg$p), size = nsamp2,
                       replace = T, prob = df_marg$p/sum(df_marg$p))
samplesalphabeta = data.frame(alpha = cA[samp_indices2[1:nsamp2]], beta = cB[samp_indices2[1:nsamp2]])
scatteralphabeta = ggplot(samplesalphabeta, aes(x=alpha, y=beta)) + geom_point()
grid.arrange(scatteralphabeta, postdensityalphabeta)

transformation

samplesalphabetatransform = samplesalphabeta
samplesalphabetatransform$alpha = log(samplesalphabeta$alpha/samplesalphabeta$beta)
samplesalphabetatransform$beta = log(samplesalphabeta$alpha+samplesalphabeta$beta)
scatteralphabetatransform = ggplot(samplesalphabetatransform, aes(x=alpha, y=beta)) + geom_point() + labs(x = "log(alpha/beta)", y = "log(alpha+beta)")
print(summary(samplesalphabetatransform$alpha))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -2.159  -1.855  -1.783  -1.782  -1.712  -1.384
print(summary(samplesalphabetatransform$beta))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.677   2.520   2.739   2.742   2.963   3.646
scatteralphabetatransform

calculate beta density at samples of p(alpha,beta|data)

df_psamp <- mapply(function(a, b, x) dbeta(x, a, b),
                   samp_A, samp_B, MoreArgs = list(x = x)) %>%
  as.data.frame() %>% cbind(x) %>% gather(ind, p, -x)

# helper function to convert ind to numeric for subsetting
indtonum <- function(x) strtoi(substring(x,2))

# Plot samples from the distribution of distributions Beta(alpha,beta),
# that is, plot Beta(alpha,beta) using posterior samples of alpha and beta
title2 <- 'Posterior samples from the distribution of distributions Beta(a,b)'
plot_psamp <- ggplot(data = subset(df_psamp, indtonum(ind) <= 20)) +
  geom_line(aes(x = x, y = p, group = ind)) +
  labs(x = expression(theta), y = '', title = title2) +
  scale_y_continuous(breaks = NULL)
plot_psamp

calculate predictive distribution for new theta

# The average of above distributions, is the predictive distribution
# for a new theta, and also the prior distribution for theta_j
df_psampmean <- spread(df_psamp, ind, p) %>% subset(select = -x) %>%
  rowMeans() %>% data.frame(x = x, p = .)

title3 <- 'Predictive distribution for a new theta and prior for theta_j'
plot_psampmean <- ggplot(data = df_psampmean) +
  geom_line(aes(x = x, y = p)) +
  labs(x = expression(theta), y = '', title = title3) +
  scale_y_continuous(breaks = NULL)
plot_psampmean

# combine the plots
grid.arrange(plot_psamp, plot_psampmean)

compare separate model and hierarchical model

# And finally compare the separate model and hierarchical model
# (using every seventh sample for clarity)

# the separate model
plot_sep7 <- ggplot(data = subset(df_sep, indtonum(ind)%%7==0)) +
  geom_line(aes(x = x, y = p, color = (ind=='V49'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Separate model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue', 'red'), guide = F) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))

# the hierarchical model
# Note that these marginal posteriors for theta_j are more narrow than
# in the separate model case, due to the borrowed information from the
# other theta_j's

# average density over samples (of a and b) for each (n,y)-pair
# at each point x
bdens2 <- function(n, y, a, b, x)
  rowMeans(mapply(dbeta, a + y, n - y + b, MoreArgs = list(x = x)))

df_hier <- mapply(bdens2, n, y, MoreArgs = list(samp_A, samp_B, x)) %>%
  as.data.frame() %>% cbind(x) %>% gather(ind, p, -x)

plot_hier7 <- ggplot(data = subset(df_hier, indtonum(ind)%%7==0)) +
  geom_line(aes(x = x, y = p, color = (ind=='V49'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Hierarchical model', color = '') +
  scale_color_manual(values = c('blue', 'red'), guide = F) +
  scale_y_continuous(breaks = NULL) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))

grid.arrange(plot_sep7, plot_hier7)

use posterior samples of theta to approximate

closed form (beta) posterior distributions of theta

# sample from posterior distribution of theta_j

nsamp <- 10000
samp_indices <- sample(length(df_marg$p), size = nsamp,
                       replace = T, prob = df_marg$p/sum(df_marg$p))

samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]
samplestheta <- matrix(0, nsamp, length(y))
for(j in 1:length(y)){
  samplestheta[, j] = sapply(1:nsamp, function(k) rbeta(1, samp_A[k]+y[j], samp_B[k]+n[j]-y[j]))
}
breakpoints <- seq(0, 1, length.out = 200)
xx <- hist(samplestheta[, 1], breaks = breakpoints, plot = FALSE)$mids
densitytheta <- matrix(0, length(breakpoints) - 1, length(y))
for(j in 1:length(y)){
  densitytheta[, j] <- hist(samplestheta[, j], breaks = breakpoints, plot = FALSE)$density
}
df_hier_samp <-densitytheta %>%
  as.data.frame() %>% cbind(xx) %>% gather(ind, p, -xx)
plot_hier7_samp <- ggplot(data = subset(df_hier_samp, indtonum(ind)%%7==0)) +
  geom_line(aes(x = xx, y = p, color = (ind=='V49'), group = ind)) +
  labs(x = expression(theta), y = '', title = 'Hierarchical model', color = '') +
  scale_color_manual(values = c('blue', 'red'), guide = F) +
  scale_y_continuous(breaks = NULL) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
grid.arrange(plot_hier7_samp, plot_hier7)

making inference on theta’s

nsamp <- 1000
samp_indices <- sample(length(df_marg$p), size = nsamp,
                       replace = T, prob = df_marg$p/sum(df_marg$p))

samp_A <- cA[samp_indices[1:nsamp]]
samp_B <- cB[samp_indices[1:nsamp]]
samplestheta <- matrix(0, nsamp, length(y))
for(j in 1:length(y)){
  samplestheta[, j] = sapply(1:nsamp, function(k) rbeta(1, samp_A[k]+y[j], samp_B[k]+n[j]-y[j]))
}
postintervals_sample <- apply(samplestheta, 2, function(x) c(median(x), c(quantile(x, c(0.025, 0.975)))))

visualize results and compare

rawest <- jitter(y / n)
plot(rawest, postintervals_sample[1, ], pch = 19, col = 'red', cex = 0.5, ylim = c(-0.01, 0.4), ylab = 'posterior intervals', xlab = 'raw estimator')
for(k in 1:length(y)){
  lines(cbind(rep(rawest[k], 2), postintervals_sample[2:3, k]))
}
lines(seq(0, 0.4, by = 0.01), seq(0, 0.4, by = 0.01), col = 'blue')
phatpool = sum(y)/sum(n)
points(phatpool, phatpool, cex = 2, col = 'green')
lines(cbind(rep(qbeta(0.5, sum(y)+1, sum(n)-sum(y)+1), 2), c(qbeta(0.025, sum(y)+1, sum(n)-sum(y)+1), qbeta(0.975, sum(y)+1, sum(n)-sum(y)+1))), col = 'yellow', lwd = 4)

Eight schools example

# Hierarchical model for SAT-example data 

y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)

# Plot data, use normpdf to plot both the y_j and sigma_j
x <- seq(-40, 60, length.out = 500)

df_sep <- mapply(function(y, s, x) dnorm(x, y, s), y, s, MoreArgs = list(x = x)) %>%
  as.data.frame() %>% setNames(LETTERS[1:8]) %>% cbind(x) %>% gather(school, p, -x)

labs1 <- c('Other Schools', 'School A')
plot_sep <- ggplot(data = df_sep) +
  geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
  labs(x = 'Treatment effect', y = '', title = 'Separate model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
plot_sep

pooled estimate

df_pool <- data.frame(x = x, p = dnorm(x, sum(y/s^2)/sum(1/s^2), sqrt(1/sum(1/s^2))))

plot_pool <- ggplot(data = df_pool) +
  geom_line(aes(x = x, y = p, color = '1')) +
  labs(x = 'Treatment effect', y = '', title = 'Pooled model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = 'red', labels = 'All schools') +
  theme(legend.background = element_blank(), legend.position = c(0.7,0.9))
plot_pool

ggplot(data = df_sep) +
  geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
  labs(x = 'Treatment effect', y = '', title = 'Separate vs pooled model', color = '') +
  scale_y_continuous(breaks = NULL) +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.8,0.9))+
  geom_line(data = df_pool, aes(x = x, y = p), color = 'green')

## conditional posteriors theta|tau, y

y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)
condmean <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  condmeanj <- (hatmu/tau^2 + y / s^2) / (1/s^2 + 1/tau^2)
  return(condmeanj)
}
condvar <- function(tau){
  sigtau <- s^2 + tau^2
  Vmu <- 1/sum(1/sigtau)
  condvarj <- 1 / (1/s^2 + 1/tau^2) + 1/tau^4 / (1/tau^2 + 1/s^2)^2 * Vmu
}
taurange <- seq(0.01, 30, length.out = 100)
condmeantau <- sapply(taurange, condmean)
condvartau <- sapply(taurange, condvar)
condmeantaudf <- t(condmeantau) %>% as.data.frame() %>% setNames(LETTERS[1:8]) %>% cbind(taurange) %>% gather(school, p, -taurange)

plot_condmeans <- ggplot(data = condmeantaudf) +
  geom_line(aes(x = taurange, y = p, color = (school=='A'), group = school)) +
  coord_cartesian(ylim = range(condmeantaudf$p) + c(-4, 5)) +
  labs(x = expression(tau), y = 'E[theta_j|tau,y)', title = '', color = '') +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.1,0.9))

condvartaudf <- t(condvartau) %>% as.data.frame() %>% setNames(LETTERS[1:8]) %>% cbind(taurange) %>% gather(school, p, -taurange)

plot_condvars <- ggplot(data = condvartaudf) +
  geom_line(aes(x = taurange, y = p, color = (school=='A'), group = school)) +
  coord_cartesian(ylim = range(condvartaudf$p) + c(-4, 5)) +
  labs(x = expression(tau), y = 'Var[theta_j|tau,y)', title = '', color = '') +
  scale_color_manual(values = c('blue','red'), labels = labs1) +
  theme(legend.background = element_blank(), legend.position = c(0.1,0.9))

grid.arrange(plot_condmeans, plot_condvars)

Marginal posterior distribution of tau

y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)
logpost_tau <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  L <- sum(log(sigtau)) / 2 + sum((y - hatmu)^2/(2 * sigtau)) - log(Vmu) / 2
  return(-L)
}
taurange <- seq(0, 30, length.out = 1000)
denstau <- sapply(taurange, logpost_tau)
denstau <- exp(denstau - max(denstau))
denstau <- denstau / (sum(denstau) * (taurange[2] - taurange[1]))
plot(taurange, denstau, xlab = expression(tau), ylab = 'marginal posterior density', main = 'marginal posterior of tau')

## sample from posterior

Nsample <- 5000
# first sample tau
y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)
logpost_tau <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  L <- sum(log(sigtau)) / 2 + sum((y - hatmu)^2/(2 * sigtau)) - log(Vmu) / 2
  return(-L)
}
taurange <- seq(0.01, 30, length.out = 1000)
denstau <- sapply(taurange, logpost_tau)
tausamples <- sample(taurange, replace = TRUE, prob = exp(denstau - max(denstau)))
# then sample mu
calculatemeanva <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  return(c(hatmu, Vmu))
}
meanvar <- sapply(tausamples, calculatemeanva)
samplesmu = apply(meanvar, 2, function(x) rnorm(1, x[1], sqrt(x[2])))
# then sample thetaj
thetasamples <- matrix(0, Nsample, length(y))
for(j in 1:length(y)){
  thetaj = (y[j] /s[j]^2 + samplesmu / tausamples^2) / (1/s[j]^2 + 1/tausamples^2)
  varj = 1/(1/s[j]^2 + 1/ tausamples^2)
  thetasamples[, j] = rnorm(Nsample, thetaj, sqrt(varj))
}

summarize posterior samples

print(apply(thetasamples, 2, function(x) quantile(x, c(0.025, 0.25, 0.5, 0.75, 0.975))))
##            [,1]      [,2]       [,3]      [,4]      [,5]      [,6]
## 2.5%  -1.825053 -4.559723 -10.660599 -5.287360 -8.277041 -8.046812
## 25%    6.291572  4.267319   2.504476  4.011129  1.874727  2.386996
## 50%   10.474975  8.066716   6.881305  7.810412  5.959016  6.766654
## 75%   15.573270 11.957768  11.018584 11.890583  9.643603 10.690779
## 97.5% 31.263886 20.016015  20.454407 20.770849 16.941934 18.926983
##            [,7]      [,8]
## 2.5%  -1.084632 -6.363821
## 25%    6.246072  4.372424
## 50%   10.156009  8.355243
## 75%   14.442019 12.774814
## 97.5% 26.133299 25.095802
par(mfrow = c(1, 2))
hist(thetasamples[, 1], main = 'School A', xlim = c(-20, 60), xlab = 'theta_j', ylab = '', freq = FALSE)
hist(apply(thetasamples, 1, max), xlim = c(-20, 60), main = 'Max', xlab = 'max theta', ylab = '', freq = FALSE)

Weakly informative priors for tau

y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)
logpost_tau_flatontau <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  L <- sum(log(sigtau)) / 2 + sum((y - hatmu)^2/(2 * sigtau)) - log(Vmu) / 2
  return(-L)
}
logpost_tau_flatonlogtau <- function(tau){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  L <- sum(log(sigtau)) / 2 + sum((y - hatmu)^2/(2 * sigtau)) - log(Vmu) / 2 + log(tau)
  return(-L)
}
logpost_tau_invGamma <- function(tau, epsilon){
  sigtau <- s^2 + tau^2
  hatmu <- sum(y/sigtau) / sum(1/sigtau)
  Vmu <- 1/sum(1/sigtau)
  L <- sum(log(sigtau)) / 2 + sum((y - hatmu)^2/(2 * sigtau)) - log(Vmu) / 2 + log(tau) + epsilon / tau^2 + (epsilon + 1) * log(tau) * 2
  return(-L)
}
taurange <- seq(0.01, 1, length.out = 1000)
dens_flatontau <- sapply(taurange, logpost_tau_flatontau)
dens_flatonllogtau <- sapply(taurange, logpost_tau_flatonlogtau)
epsilon <- 1
dens_invgamma <- sapply(taurange, function(x) logpost_tau_invGamma(x, epsilon))
epsilon <- 0.1
dens_invgamma1 <- sapply(taurange, function(x) logpost_tau_invGamma(x, epsilon))
epsilon <- 0.01
dens_invgamma2 <- sapply(taurange, function(x) logpost_tau_invGamma(x, epsilon))
epsilon <- 0.001
dens_invgamma3 <- sapply(taurange, function(x) logpost_tau_invGamma(x, epsilon))
dens_flatontau <- exp(dens_flatontau - max(dens_flatontau))
dens_flatonllogtau <- exp(dens_flatonllogtau - max(dens_flatonllogtau))
dens_invgamma <- exp(dens_invgamma - max(dens_invgamma))
dens_invgamma1 <- exp(dens_invgamma1 - max(dens_invgamma1))
dens_invgamma2 <- exp(dens_invgamma2 - max(dens_invgamma2))
dens_invgamma3 <- exp(dens_invgamma3 - max(dens_invgamma3))
dens_flatontau <- dens_flatontau / sum(dens_flatontau)
dens_flatonllogtau <- dens_flatonllogtau / sum(dens_flatonllogtau)
dens_invgamma <- dens_invgamma / sum(dens_invgamma)
dens_invgamma1 <- dens_invgamma1 / sum(dens_invgamma1)
dens_invgamma2 <- dens_invgamma2 / sum(dens_invgamma2)
dens_invgamma3 <- dens_invgamma3 / sum(dens_invgamma3)
plot(taurange, dens_flatontau, col = 'blue', ylim = range(c(dens_flatonllogtau, dens_flatontau, dens_invgamma, dens_invgamma1, dens_invgamma2, dens_invgamma3)), type = 'l', ylab = '', main = 'marginal posterior of tau')
lines(taurange, dens_flatonllogtau, col = 'red')
lines(taurange, dens_invgamma, col = 'green')
lines(taurange, dens_invgamma1, col = 'black')
lines(taurange, dens_invgamma2, col = 'magenta')
lines(taurange, dens_invgamma3, col = 'yellow')